<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final//EN">
<!--Converted with LaTeX2HTML 98.2 beta6 (August 14th, 1998)
original version by:  Nikos Drakos, CBLU, University of Leeds
* revised and updated by:  Marcus Hennecke, Ross Moore, Herb Swan
* with significant contributions from:
  Jens Lippmann, Marek Rouchal, Martin Wilck and others -->
<HTML>
<HEAD>
<TITLE>General Linear Model Problem</TITLE>
<META NAME="description" CONTENT="General Linear Model Problem">
<META NAME="keywords" CONTENT="lug_l2h">
<META NAME="resource-type" CONTENT="document">
<META NAME="distribution" CONTENT="global">
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<LINK REL="STYLESHEET" HREF="lug_l2h.css">
<LINK REL="previous" HREF="node85.html">
<LINK REL="up" HREF="node84.html">
<LINK REL="next" HREF="node88.html">
</HEAD>
<BODY >
<!--Navigation Panel-->
<A NAME="tex2html5410"
 HREF="node88.html">
<IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next"
 SRC="next_motif.gif"></A> 
<A NAME="tex2html5404"
 HREF="node84.html">
<IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up"
 SRC="up_motif.gif"></A> 
<A NAME="tex2html5400"
 HREF="node86.html">
<IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous"
 SRC="previous_motif.gif"></A> 
<A NAME="tex2html5406"
 HREF="node1.html">
<IMG WIDTH="65" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="contents"
 SRC="contents_motif.gif"></A> 
<A NAME="tex2html5408"
 HREF="node152.html">
<IMG WIDTH="43" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="index"
 SRC="index_motif.gif"></A> 
<BR>
<B> Next:</B> <A NAME="tex2html5411"
 HREF="node88.html">Further Details: Error Bounds</A>
<B> Up:</B> <A NAME="tex2html5405"
 HREF="node84.html">Error Bounds for Generalized</A>
<B> Previous:</B> <A NAME="tex2html5401"
 HREF="node86.html">Further Details: Error Bounds</A>
 &nbsp <B>  <A NAME="tex2html5407"
 HREF="node1.html">Contents</A></B> 
 &nbsp <B>  <A NAME="tex2html5409"
 HREF="node152.html">Index</A></B> 
<BR>
<BR>
<!--End of Navigation Panel-->

<H2><A NAME="SECTION03462000000000000000">
General Linear Model Problem</A>
</H2>

<P>
The general linear model (GLM) problem is
<BR><P></P>
<DIV ALIGN="CENTER">

<!-- MATH
 \begin{displaymath}
\min_{x,y} \|y\| \quad \mbox{subject to} \quad d = Ax + By
\end{displaymath}
 -->


<IMG
 WIDTH="269" HEIGHT="40" BORDER="0"
 SRC="img491.gif"
 ALT="\begin{displaymath}\min_{x,y} \Vert y\Vert \quad \mbox{subject to} \quad d = Ax + By \end{displaymath}">
</DIV>
<BR CLEAR="ALL">
<P></P>
where <B><I>A</I></B> is an <B><I>n</I></B>-by-<B><I>m</I></B> matrix, <B><I>B</I></B> is an <B><I>n</I></B>-by-<B><I>p</I></B> matrix, and <B><I>d</I></B>
is a given <B><I>n</I></B>-vector, with 
<!-- MATH
 $m \leq n \leq m+p$
 -->
<IMG
 WIDTH="122" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
 SRC="img21.gif"
 ALT="$m \leq n \leq m+p$">.
<BR>

<P>
The GLM problem is solved by the driver routine xGGGLM
(see section <A HREF="node84.html#sec_lseglm_drivers">4.6</A>).
Let <IMG
 WIDTH="14" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
 SRC="img467.gif"
 ALT="$\widehat {x}$">
and <IMG
 WIDTH="14" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
 SRC="img492.gif"
 ALT="$\widehat {y}$">
be the computed values of <B><I>x</I></B> and <B><I>y</I></B>, respectively.
The approximate error bounds
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG
 WIDTH="441" HEIGHT="111" BORDER="0"
 SRC="img493.gif"
 ALT="\begin{eqnarray*}
% latex2html id marker 11027\frac{ \Vert x - \widehat {x} \...
...hfil ...">
</DIV><P></P>
<BR CLEAR="ALL">
can be computed by the following code fragment

<P>
<PRE>
      EPSMCH = SLAMCH( 'E' )
*     Compute the 2-norm of the left hand side D
      DNORM = SNRM2( N, D, 1 )
*     Solve the generalized linear model problem
      CALL SGGGLM( N, M, P, A, LDA, B, LDB, D, Xc, Yc, WORK,
     $            LWORK, IWORK, INFO )
*     Compute the F-norm of A and B
      ANORM = SLANTR( 'F', 'U', 'N', M, M, A, LDA, WORK( M+NP+1 ) )
      BNORM = SLANTR( 'F', 'U', 'N', N, P, B( 1, MAX( 1, P-N+1 ) ),
     $                LDB, WORK( M+NP+1 ) )
*     Compute the 2-norm of Xc
      XNORM = SNRM2( M, Xc, 1 )
*     Condition estimation
      IF( N.EQ.M ) THEN
         PBPSNM = ZERO
         TNORM = SLANTR( '1', 'U', 'N', N, N, A, LDA, WORK( M+NP+M+1 ) )
         CALL STRCON( '1', 'U', 'N', N, A, LDA, RCOND, WORK( M+NP+M+1 ),
     $                IWORK, INFO )
         ABPSNM = ONE / (RCOND * TNORM )
      ELSE
*        Compute norm of (PB)^+
         TNORM = SLANTR( '1', 'U', 'N', N-M, N-M, B( M+1, P-N+M+1 ), LDB,
     $                   WORK( M+NP+1 ) )
         CALL STRCON( '1', 'U', 'N', N-M, B( M+1, P-N+M+1 ), LDB, RCOND,
     $                WORK( M+NP+1 ), IWORK, INFO )
         PBPSNM = ONE / (RCOND * TNORM )
*        Estimate norm of A^+_B
         KASE = 0
         CALL SLACON( N, WORK( M+NP+1 ), WORK( M+NP+N+1 ), IWORK, EST, KASE )
   30    CONTINUE
            CALL STRSV( 'Upper', 'No transpose', 'Non unit', N-M,
     $                  B( M+1, P-N+M+1 ), LDB, WORK( M+NP+N+M+1 ), 1 )
            CALL SGEMV( 'No transpose', M, N-M, -ONE, B( 1, P-N+M+1 ),
     $                  LDB, WORK( M+NP+N+M+1 ), 1, ONE,
     $                  WORK( M+NP+N+1 ), 1 )
            CALL STRSV( 'Upper', 'No transpose', 'Non unit', M, A, LDA,
     $                  WORK( M+NP+N+1 ), 1 )
            DO I = 1, P
               WORK( M+NP+I ) = WORK( M+NP+N+I )
            END DO
            CALL SLACON( M, WORK( M+NP+N+1 ), WORK( M+NP+1 ), IWORK, EST, KASE )
            IF( KASE.EQ.0 ) GOTO 40
            CALL STRSV( 'Upper', 'Transpose', 'Non unit', M, A, LDA,
     $                  WORK( M+NP+1 ), 1 )
            CALL SGEMV( 'Transpose', M, N-M, -ONE, B( 1, P-N+M+1 ), LDB,
     $                  WORK( M+NP+1 ), 1, ZERO, WORK( M+NP+M+1 ), 1 )
            CALL STRSV( 'Upper', 'Transpose', 'Non unit', N-M,
     $                  B( M+1, P-N+M+1 ), LDB, WORK( M+NP+M+1 ), 1 )
            DO I = 1, N
               WORK( M+NP+N+I ) = WORK( M+NP+I )
            END DO
            CALL SLACON( N, WORK( M+NP+1 ), WORK( M+NP+N+1 ), IWORK, EST, KASE )
            IF( KASE.EQ.0 ) GOTO 40
         GOTO 30
   40    CONTINUE
         ABPSNM = EST
      END IF
*     Estimate norm of (A^+_B)*B
      IF( P+M.EQ.N ) THEN
         EST = ZERO
      ELSE
         KASE = 0
         CALL SLACON( P-N+M, WORK( M+NP+1 ), WORK( M+NP+M+1 ), IWORK, EST, KASE )
   50    CONTINUE
*
            IF( P.GE.N ) THEN
               CALL STRMV( 'Upper', 'No trans', 'Non Unit', M,
     $                     B( 1, P-N+1 ), LDB, WORK( M+NP+M+P-N+1 ), 1 )
               DO I = 1, M
                  WORK( M+NP+I ) = WORK( M+NP+M+P-N+I )
               END DO
            ELSE
               CALL SGEMV( 'No transpose', N-P, P-N+M, ONE, B, LDB,
     $                  WORK( M+NP+M+1 ), 1, ZERO, WORK( M+NP+1 ), 1 )
               CALL STRMV( 'Upper', 'No trans', 'Non Unit', P-N+M,
     $                     B( N-P+1, 1 ), LDB, WORK( M+NP+M+1 ), 1 )
               DO I = N-P+1, M
                  WORK( M+NP+I ) = WORK( M+NP+M-N+P+I )
               END DO
            END IF
            CALL STRSV( 'Upper', 'No transpose', 'Non unit', M, A, LDA,
     $                  WORK( M+NP+1 ), 1 )
            CALL SLACON( M, WORK( M+NP+M+1 ), WORK( M+NP+1 ), IWORK, EST, KASE )
*
            IF( KASE.EQ.0 ) GOTO 60
*
            CALL STRSV( 'Upper', 'Transpose', 'Non unit', M, A, LDA,
     $                  WORK( M+NP+1 ), 1 )
            IF( P.GE.N ) THEN
               CALL STRMV( 'Upper', 'Trans', 'Non Unit', M,
     $                     B( 1, P-N+1 ), LDB, WORK( M+NP+1 ), 1 )
               DO I = 1, M
                  WORK( M+NP+M+P-N+I ) = WORK( M+NP+I )
               END DO
               DO I = 1, P-N
                  WORK( M+NP+M+I ) = ZERO
               END DO
            ELSE
               CALL STRMV( 'Upper', 'Trans', 'Non Unit', P-N+M,
     $                     B( N-P+1, 1 ), LDB, WORK( M+NP+N-P+1 ), 1 )
               DO I = 1, P-N+M
                  WORK( M+NP+M+I ) = WORK( M+NP+N-P+I )
               END DO
               CALL SGEMV( 'Transpose', N-P, P-N+M, ONE, B, LDB,
     $                  WORK( M+NP+1 ), 1, ONE, WORK( M+NP+M+1 ), 1 )
            END IF
            CALL SLACON( P-N+M, WORK( M+NP+1 ), WORK( M+NP+M+1 ), IWORK, EST,
     $                   KASE )
*
            IF( KASE.EQ.0 ) GOTO 60
         GOTO 50
   60    CONTINUE
      END IF
      ABPSBN = EST
*     Get condition numbers and approximate error bounds
      CNDAB = ANORM*ABPSNM
      CNDBA = BNORM*PBPSNM
      IF( PBPSNM.EQ.0.0E+0 ) THEN
*        Then A is square and nonsingular
         XERRBD = EPSMCH*( CNDAB*( ONE+DNORM/(ANORM*XNORM) ) )
         YERRBD = 0.0E+0
      ELSE
         XERRBD = EPSMCH*( CNDAB*( ONE+DNORM/(ANORM*XNORM) ) +
     $                2.0E0*CNDAB*CNDBA*CNDBA*DNORM/(ANORM*XNORM) +
     $                ABPSBN*ABPSBN*PBPSNM*PBPSNM*ANORM*DNORM/XNORM )
         YERRBD = EPSMCH*( ABPSBN*ANORM*PBPSNM*PBPSNM +
     $                PBPSNM*(ANORM*XNORM/DNORM + 2.0E0*CNDBA*CNDBA +
     $                ONE) + CNDBA*PBPSNM )
      END IF
</PRE>

<P>
For example, if 
<!-- MATH
 ${\tt SLAMCH('E')} = 2^{-24} = 5.961 \cdot 10^{-8}$
 -->
<IMG
 WIDTH="259" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
 SRC="img397.gif"
 ALT="${\tt SLAMCH('E')} = 2^{-24} = 5.961 \cdot 10^{-8}$">,
<BR><P></P>
<DIV ALIGN="CENTER">

<!-- MATH
 \begin{displaymath}
A = \left( \begin{tabular}{rrrr}
                  1 &  2 &  1 & 4 \\
                  1 &  3 &  2 & 1 \\
                 -1 & -2 & -1 & 1 \\
                 -1 &  2 & -1 & 5 \\
                  1 &  0 &  0 & 1
               \end{tabular} \right), \quad
    B = \left( \begin{tabular}{rrr}
                  1 &  2 &  2 \\
                 -1 &  1 & -2\\
                  3 &  1 &  6 \\
                  1 & -1 &  2 \\
                  2 & -2 &  4
               \end{tabular} \right)
    \quad \mbox{and} \quad
    d = \left( \begin{array}{r}
                   1 \\
                   1 \\
                   1 \\
                   1 \\
                   1
               \end{array} \right).
\end{displaymath}
 -->


<IMG
 WIDTH="544" HEIGHT="115" BORDER="0"
 SRC="img494.gif"
 ALT="\begin{displaymath}A = \left( \begin{tabular}{rrrr}
1 &amp; 2 &amp; 1 &amp; 4 \\
1 &amp; 3 &amp; ...
...array}{r}
1 \\
1 \\
1 \\
1 \\
1
\end{array} \right). \end{displaymath}">
</DIV>
<BR CLEAR="ALL">
<P></P>
Then (to 7 decimal places)
<BR><P></P>
<DIV ALIGN="CENTER">

<!-- MATH
 \begin{displaymath}
\widehat {x} = \left( \begin{array}{r}
                            -0.5466667  \\
                             0.3200001  \\
                             0.7200000  \\
                            -0.0533334
                        \end{array} \right)
    \quad \mbox{and} \quad
    \widehat {y} = \left( \begin{array}{r}
                         0.1333334 \\
                        -0.1333334 \\
                         0.2666667
                    \end{array} \right).
\end{displaymath}
 -->


<IMG
 WIDTH="402" HEIGHT="93" BORDER="0"
 SRC="img495.gif"
 ALT="\begin{displaymath}\widehat {x} = \left( \begin{array}{r}
-0.5466667 \\
0.320...
... 0.1333334 \\
-0.1333334 \\
0.2666667
\end{array} \right).\end{displaymath}">
</DIV>
<BR CLEAR="ALL">
<P></P>
The computed error bounds

<!-- MATH
 ${\tt XERRBD} = 1.2\cdot10^{-5}$
 -->
<IMG
 WIDTH="153" HEIGHT="19" ALIGN="BOTTOM" BORDER="0"
 SRC="img496.gif"
 ALT="${\tt XERRBD} = 1.2\cdot10^{-5}$">
and 
<!-- MATH
 ${\tt YERRBD} = 6.4\cdot10^{-7}$
 -->
<IMG
 WIDTH="153" HEIGHT="19" ALIGN="BOTTOM" BORDER="0"
 SRC="img497.gif"
 ALT="${\tt YERRBD} = 6.4\cdot10^{-7}$">,
where 
<!-- MATH
 ${\tt CNDAB} = 26.97$
 -->
<IMG
 WIDTH="113" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
 SRC="img498.gif"
 ALT="${\tt CNDAB} = 26.97$">,

<!-- MATH
 ${\tt CNDBA} = 1.78$
 -->
<IMG
 WIDTH="104" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
 SRC="img499.gif"
 ALT="${\tt CNDBA} = 1.78$">,
The true errors in <B><I>x</I></B> and <B><I>y</I></B> are 
<!-- MATH
 $1.2\cdot10^{-7}$
 -->
<IMG
 WIDTH="75" HEIGHT="19" ALIGN="BOTTOM" BORDER="0"
 SRC="img500.gif"
 ALT="$1.2\cdot10^{-7}$">
and

<!-- MATH
 $1.3\cdot10^{-7}$
 -->
<IMG
 WIDTH="75" HEIGHT="19" ALIGN="BOTTOM" BORDER="0"
 SRC="img501.gif"
 ALT="$1.3\cdot10^{-7}$">,
respectively.  Note that the exact solutions are

<!-- MATH
 $x = \frac{1}{75}[-41,24,54,-4]^T$
 -->
<IMG
 WIDTH="187" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
 SRC="img502.gif"
 ALT="$x = \frac{1}{75}[-41,24,54,-4]^T$">
and

<!-- MATH
 $y = \frac{1}{15}[2,-2,4]^T$
 -->
<IMG
 WIDTH="130" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
 SRC="img503.gif"
 ALT="$y = \frac{1}{15}[2,-2,4]^T$">.

<P>
<BR><HR>
<!--Table of Child-Links-->
<A NAME="CHILD_LINKS"></A>

<UL>
<LI><A NAME="tex2html5412"
 HREF="node88.html">Further Details:  Error Bounds for General Linear Model Problems</A>
</UL>
<!--End of Table of Child-Links-->
<BR><HR>
<ADDRESS>
<I>Susan Blackford</I>
<BR><I>1999-10-01</I>
</ADDRESS>
</BODY>
</HTML>
